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Abstract 

The formation and growth of hehum bubbles due to self-irradiation in plutonium has 
been modelled by a discrete kinetic equations for the number densities of bubbles 
having k atoms. Analysis of these equations shows that the bubble size distribution 
function can be approximated by a composite of: (i) the solution of partial differen- 
tial equations describing the continuum limit of the theory but corrected to take into 
account the effects of discreteness, and (ii) a local expansion about the advancing 
leading edge of the distribution function in size space. Both approximations con- 
tribute to the memory term in a close integrodifferential equation for the monomer 
concentration of single helium atoms. The present boundary layer theory for dis- 
crete equations is compared to the numerical solution of the full kinetic model and 
to previous approximation of Schaldach and Wolfer involving a truncated system of 
moment equations. 
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1 Introduction 



There are simple kinetic models of irreversible aggregation, in which a cluster 
with k monomers grows by absorbing one monomer but it cannot decrease 
in size by shedding part of its mass. An interesting example is the forma- 
tion and growth of helium bubbles in plutonium alloys as a consequence of 
alpha decay due to self-irradiation [1,2]. As an alloy ages, there is an initial 
transient stage during which self-irradiation produces dislocation loops that 
tend to saturate within approximately two years. The alpha particles created 
during irradiation become helium atoms. These atoms come to rest at unfilled 
vacancies generated during their slowing-down process, before they are cap- 
tured at existing helium bubbles. A helium atom diffuses through the lattice 
until it finds another helium atom thereby forming a stable dimer or until it 
finds a helium bubble (a stable cluster with k atoms or, in short, a A;-cluster), 
which absorbes it. Helium bubbles are attached to lattice defects, do not move 
and do not shed helium atoms because the binding energies of helium to any 
cluster are extremely high [2]. The following kinetic model based on these 
observations has been proposed by Schaldach and Wolfer [1]: 



Pk = 4:7tD Pi Gk-iPk-i - 47rL' pi akpk, k > 3, (1) 
P2 = SttD plai - AttD pia2P2, (2) 

oo * 

Pi + T.kpk= g{t') dt'. (3) 

k=2 i 

Here pk = dpk/dt, pk is the number density of k clusters having effective radii 
Ofc (when the center of a monomer comes within distance of the cluster 
center, it is absorbed), pi is the number of monomers per unit volume, D is 
the diffusion coefficient and g(t) is the number of monomers created per unit 
volume and per unit time. Eq. (3) means that the total number of monomers 
per unit volume, whether they are in solution or forming part of a A; cluster, 
should equal the time integral of g(i). In Equations (1) and (2), k clusters 
grow by adding one monomer with a rate AirDpiak-i (for k > 2) to a, k — 1 
cluster, and they do not decay. The mean absorption rate of a monomer by 
an immobile k cluster (with A; > 1) is AnDpiak, whereas the rate of creation 
of an immobile dimer by the collision of two mobile monomers is twice this 
quantity, SnDpiai. The absorption rate can be calculated immediately from 
the concentration field outside a k cluster (i.e., the number of monomers per 
unit volume in f > a^), p{r,i) assuming that the mean-field approximation 
holds. The process of diffusion-limited adiabatic growth of clusters is random: 
diffusion of monomers to a given cluster is randomly affected by the spatial 
distribution of the other clusters in its vicinity. If alpk -C 1 (diluted clusters), 
we may assume that the concentration field outside a cluster goes to the 
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average monomer density in solution, pi{i), which ignores fluctuations in the 
local environment and constitutes the mean-field approximation. Then p{f, t) 
solves the Laplace equation outside f = a^, with p{ak,i) = 0, and it takes on 
the value pi(t) as f — > oo. Thus 

p{f,i)^p^{i) (l-y), 
for r > a^- The mean absorption rate of monomers by the k cluster is 



inalD-^iak, i) = ATrDp^ak- (4) 

The simplest model for the efi'ective radius of a A; cluster is based on packing 
of non-overlapping particles: 



(5) 



The mean-field model (1) - (5) has obvious limitations: the kinetic equations 
are deterministic due to the mean-field approximation while the creation of 
helium atoms and diffusion of monomers are random processes, and helium 
bubbles have limitations to their growth beyond a certain size. These cff'ccts 
will not be addressed in the present paper. It is interesting to observe that a 
related kinetic system was proposed and solved in 1914 by McKendrick as a 
model of leucocyte phagocytosis [3]. In McKendrick's model, pk is the density 
of leucocytes which have ingested k bacteria, and its rate equation is (1) for 
k > 0, with a known function of time pi{t) > and p_i = 0. McKendrick's 
solution method involved solving the equation for po in terms of / pi dt and 
solving recursively all other equations for p^ as functions of po- His method 
cannot be used to solve the system (1) - (3), but an useful closure of this 
infinite system to only three differential equations was introduced in [1], and 
compared to experiments, [1,2]. 

In this paper, we shall study the solution of the system (1) - (3) starting from 
an initial condition corresponding to the absence of helium bubbles, i.e., 

Pfc(O) = 0, forA; > 1. (6) 

We shall consider the case of a constant production rate of helium atoms, 
g{t) = gi. Once pi(t) as been found, the total density of bubbles can be 
obtained by integrating the equation: 



d °° 

-J2pk^8nDa,pl (7) 

k=2 
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Eq. (7) is immediately obtained by adding (2) and (1) (for all A; > 3). 

We have found that, after a short transient, the size distribution function 
becomes a smooth function of k and it is a functional of the monomer con- 
centration. In turn, the monomer concentration satisfies an integrodifferential 
equation. The size distribution function can be approximated by matched 
asymptotic expansions. Except for large k, its form is close to the solution 
of the hyperbolic equation resulting from replacing derivatives instead of dif- 
ferences in (1), but including corrections due to discreteness. This outer ap- 
proximation breaks down at large k: it has an unrealistic large peak at a 
maximum vahic of k after which the distribution function abruptly falls to 
zero. If our governing equations were partial differential equations instead of 
differential-difference equations, we would say that a boundary layer should 
be inserted to remedy this imperfection. However, the concept of boundary 
layer is not straightforward for discrete equations. How do we insert a bound- 
ary layer approximation to difference equations? The answer to this question 
lies in the wave front expansion of the Beckcr-Doring equations some of us 
introduced in Ref . [4] . This theory yields one equation for the leading edge of 
the size distribution function and one equation for the distribution function 
itself near its edge. A solution of the latter written in similarity variables is 
then matched to the outer approximation and its contribution to the integral 
terms of the equation for the monomer concentration is calculated. These two 
steps were not needed for the nucleation calculation in Ref. [4] because there 
the outer approximation was a constant, size-independent profile and the in- 
tegral condition was identically satisfied during the nucleation transient. The 
present overall theory including outer and boundary layer approximations de- 
scribes quite well the observed numerical simulations of the discrete system of 
equations. 

The rest of the paper is as follows. The nondimensional equations of the model 
are derived in Section 2. The outer approximation to the size distribution func- 
tion and the monomer concentration is described in Section 3. The boundary 
layer analysis is given in Section 4. Sec. 5 compares the present theory to a 
simple system of three differential equations for the helium density, the bub- 
ble density and the monomer concentration obtained by closing the moment 
equations for the size distribution according to a simple ansatz introduced 
in [1]. The last section summarizes our conclusions, and technical details are 
relegated to Appendices. 



2 Nondimensionalization of the kinetic equations 

Let [t], [k], [c], [p] be typical units of time, cluster size, monomer concentration 
and cluster density, respectively. Notice that we use [c] instead of the more 
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cumbersome notation [pi] in order to avoid confusion between the units of 
monomer concentration and cluster density. Equations (1), (7) and (3) provide 
the following dominant balances between these scaling units: 



j^^A^D[c]a,-^, (8) 

lP^^4nD[cra,, (9) 

mp]=g\i]. (10) 

To derive these equations, we have assumed that there is a continuum limit 
with k 1 and that the number density of particles in clusters, ^,^=2^ Pki 
is much larger than the monomer concentration pi. Note that we have three 
relations (8) - (10) for four unknowns [t], [k\, [c], [p\. In the absence of another 
relation (for example, a different creation rate g{t) — G{t/T), which fixes the 
time scale [t] = r), we can express three of these scaling units in terms of the 
fourth. We find 

ra=vS5- H = /^W-''^- W = /=|=^W-" (11) 

We shall now define nondimensional units as 



Pi _ 1/2 
C—YT = 1^ Pi 



l^nDai 



[c]- ^^V 9 



Pk _ 5/6 "^T^Dai 

'''^W]^ (12) 

in which we have set [k] = n (any positive number). Equations (1) - (3) and 
(7) become 



±t^cK'/'[{k-i)'/'n.,-k'/'nl k>3, (13) 

^^2Kc'-K'/'2'^'cr2, (14) 

oo 

K-^/^C+ K-^Y.k'T'k^t, (15) 
fc=2 

|f:r, = 2«c^. (16) 

k=2 
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Our nondimcnsional system of kinetic equations comprises (13) - (15) with 
K = 1 and initial conditions r;^.(0) = 0, c(0) = 0. (16) is a consequence of the 
previous equations or it can be used instead of (14). Defining an adaptive time 
s — /q c{t') dt', we can rewrite our system in the more convenient form: 



±t ^ - ly/'r,., - k'/'n, k>3, (17) 
as 

^ = 20-2^/^2, (18) 
ds 

oo 

c+J2krk = t, (19) 

k=2 

ds , , 

P"' 

with initial conditions c(0) = rfe(O) = (A; > 2) and t(0) ~ 0. In order to 
integrate this system of equations, it is convenient to time differentiate (19) 
and use (20). The result is 



dc 

c^ + 4c2 + cMv3 = l, (21) 
in which we have defined the moments of the size distribution function rfe(s): 



= E k^^k. (22) 

k=2 

Mq and Ml are the number densities of bubbles and of helium, respectively. 
They satisfy: 



(.3) 

dMi , , 

^ = 4c + Mi/3, (24) 

with zero initial conditions. 



The kinetic equations describing formation of helium bubbles are therefore 
(17), (18), (20), (21) and (22) with zero initial conditions. These equations 
yield the monomer concentration c and the size distribution function r^. The 
number densities of bubbles and of helium are given by (23) and (24). 
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3 Outer solution and relation to the continuum limit equations 



3.1 Continuum limit and similarity solution 

A first attempt at approximating the model equations consists of taking the 
continuum hmit of (17) - (20). In fact, assume that 

n{s) = r{k,s), (25) 
and Taylor expand r in (17) up to first order terms. The result is 

J + |(^-'V,^0, (26) 

oo 

jkrdk^t. (27) 



We have ignored the monomer concentration in (27). Integrating (26) over 
A; > 0, we obtain 

d 7 

— rdk = \im(k^^^r), 

ds J k-^O^ " 



and (23) implies the following signaling condition: 

lim(/t^/V) = 2c. (28) 
fe— >o 

The method of characteristics provides the following solution to (26) and (28) 
with initial condition c(0) = 0: 

r{k, s) - 2k-^'^ c (s - a{k)) 6 {s - a{k)) , (29) 

aik) = ^ k'/', (30) 

in which ^^(a;) = 1 if a; > and 9{x) = if x < is the Hcaviside unit step 
function. After a change of variable, the integral condition (27) becomes 



2 



2\3/2 

3 



s 

J (s - s'f/^ c{s') ds' ^ t, (31) 





or, equivalently. 
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s 

c{s) j{s-s'f'^c{s')ds' = 1. (32) 





The problem (26) - (28) and (31) has a similarity solution whose role will be 
discussed later. In fact, note that (26) - (28) are invariant under the scahng 

transformation r — ^ K~^/^r, c k~^I'^c. s k'^^^s, t — » nJ^^t^ k — > Kk, which is 
suggested by the nondimensionalization (12). Therefore the combinations x = 
k and p = s^^^r are invariant under the previous scaling transformation. 
In terms of p = p(x), (26) becomes 



dp_ p 'ix'"-i 

ds 3x Ix^'^ ' 

p = s^l^r, x = ks-^/\ (34) 
The solution of (33) is 



i?0S-5/4^-l/3 

^(x, s) = : -Tji- (35) 



(i-ix 



2/3^ 



3/4- 
+ 



Here Ro is a positive constant and f{x)+ — f{x) 6{x) for any function f{x). 
This similarity solution has integrable singularities at x = and at x = Xo = 
(2/3)^''^ (corresponding to the position of the local maximum of for large 
sizes, /cm(s), which coincides with the maximum cluster size in this case). Now 
the signaling condition (28) yields the monomer concentration: 



c{s)^^Ros-^l\ (36) 
and (20) together with s(0) = yield the relation between s and the time: 



'7 x4/7 

s^[-Rot) . (37) 



The constant Rq is found by inserting (36) in (32), which yields i?^^2/3B(3/2, 1/4) 
1, from which. 



Ro = Y{TJ^ ~ 0.837042. (38) 
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Fig. 1. Nondimensional size distribution function rfc(t) evaluated by solving the full 
model system of discrete equations (solid line) and the similarity solution (dashed 
line) at the nondimensional times 100, 200, 300, 400 and 500. 

The self-similar size distribution function has a singularity at its maximum 
size X = Xo- Fig- 1 shows that this solution is a relatively poor approximation 
to the numerical solution of the full model. To improve it, we can take into 
account the effects of discreteness by replacing 1 + 3s~^ — 3x^''^/2 instead of 
1 — 3x^/^/2 in (35). Then the singularity of the self-similar size distribution 
function moves closer to the local maximum of the exact numerical solution 
and the approximation improves, as shown in Figure 2. We observe that the 
singularity of the self-similar size distribution function occurs before the nu- 
merical size distribution function reaches its local maximum, and this effect 
becomes more noticeable as time increases. 



3.2 Outer solution for the discrete problem 

The similarity solution of the continuum equations is a poor approximation to 
the numerical solution of the discrete model. We need to correct it by including 
the effects of discreteness and by approximating better the equation for the 
monomer concentration. The corrections due to discreteness can be found by 
first solving the linear equations (17) exactly. Laplace transforming (17) and 
using rfc(O) = 0, we obtain 




(39) 
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20 40 60 80 100 120 140 160 180 
k 



Fig. 2. Same as Fig. 1 but now the dashed hne is the similarity solution corrected 
by replacing 1 + Ss'^ - 'i^l'^ jl instead of 1 - 3x^/^2 in (35). 

for k > 2, from which it follows: 



k 



i=2 



[s) = 2k-^/^ J Rk{s-s') c{s')ds'. 



(40) 
(41) 

(42) 



These equations give the exact size distribution function in terms of c(s). 
Inserting (42) in the definition of M1/3 and substituting the result in (21), we 
obtain 



c— + 4c2 + 2c [J2Rkis-s')]c{s')ds' = 1. (43) 

^==2 

This integrodifferential equation for c should be solved using the initial con- 
dition c(0) = 0. Unfortunately, it is not too useful because the integral kernel 
is written as an infinite series of the inverse Laplace transform of (41), which 
is rather unwieldly. 

To find an approximate form of Eq. (43), we assume that Rk{s) is peaked 
about its mean value: 
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In Appendix A, we show that 



aik)^h'/'-3 + ^^ {k^^). (45) 

Wc also show that the relative width of the peak in Rk{s), a/a ~ 2A;^^/^/v^ 
tends to as /c — s> oo. If we assume that c(s) is a smooth function (and this is 
not always the case), then we may approximate the kernel in the convolution 
integral (42) by Rk{s) ~ 5 (s — a(/c)), so that the size distribution function is 
given by (29) with a{k) replaced by (45). Note that the first term in the right 
hand side of (45) coincides exactly with the result (30) obtained by solving 
the continuum equations (26) and (28). Equation (43) becomes 



dr 

c^ + 4c' + 2c c{s-a{k))^l. (46) 

k=2,a{k)<s 

For large s, the position of the local maximum of the cluster size distribution 
is such that a{kM) ~ 3A;^'^/2 — 3 = s. Then we can approximate the sum in 
(46) by an integral over k from k = Q to k = kuis). Changing variables in the 
integral, (46) becomes 



c^ + 4c2 + 2^/^c f (s - s' + Sy^^ c(s') ds' ^ 1. (47) 
ds \ 3 J 



We would like to compare now the solution of (29), (45), (47) and (20) with 
zero initial conditions to the numerical solution of the exact discrete model. 
To this end, it is better to write all unknowns as functions of the time t. Our 
local theory is then 



rk{t) = 2k-^/^ c(s(t) - a{k)) e{s{t) - a{k)), 

^ + 4c^ + 2 J|c / [s{t) - s{t') + 3] V2 [c{s{t'))r dt' = 1, (48) 



ds 

dt^""' 



c(0) = s(0) = 0. 



Figure 3 shows a comparison between the monomer concentration evaluated 
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Fig. 3. (a) Monomer concentration c{t) evaluated using: (i) the numerical solution 
of the discrete equations of the model (solid line) , (ii) the local theory (48) (dashed 
line), (iii) Schaldach and Wolfer's moment equations (dot-dashed line), and (iv) the 
self-similar solution having a vertical asymptote at t = (thin solid line), (b) Same 
as in (a) for a larger range of times. All variables are written in dimensionless units. 

by solving the full discrete model equations, the local theory (48), the similar- 
ity solution and the three moment equations used by Wolfer and coworkers (cf. 
Section 5 and [1]). In Fig. 3(a), we observe that Wolfer et al's approximation 
is better for short times, but that the local theory given by (48) provides the 
best approximation as time goes to infinity, cf. Fig. 3(b). Figure 4 shows the 
size distribution function calculated at different times by using the local the- 
ory (48) (dashed lines) and the numerical solution of the full discrete model. 
The local theory yields higher cluster densities than the exact solution (which 
is also the case with the self-similar solution in Fig. 2), a much higher maxi- 
mum density but it predicts a maximum size which is very close to the local 
maximum of the real distribution function at large sizes. 



4 Leading edge of the size distribution function 



The previous local description of the size distribution function differs sub- 
stantially from the numerical solution of the model equations for large sizes. 
However, the maximum of the numerical coincides with the peak of the 
approximate ai k = kM{s), which is also its leading edge. To improve our 
asymptotic theory, we should insert a moving boundary layer there. How? We 
shall use the wave front theory which we introduced in our previous work on 
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Fig. 4. Dimensionless size distribution function rfc(t) calculated using the numerical 
solution of the full model discrete equations (solid line) and the local theory (48) 
(dashed line) at the nondimensional times 100, 200, 300, 400 and 500. 

homogeneous nucleation [4]. Firstly, let us rewrite (17) and (21) as 



^ = k'/'ia,^^-a,), k>3, 
as 

dc °° 
c— + 4c2 + c ^afc = 1. 

"'^ k=2 



(49) 
(50) 
(51) 



Secondly, let us use local coordinates about the inflection point of the wave 
front leading edge in cxfc, k = K{s), for large k: 



ak = S{X,s), X = k-K{s), 1<^X<^K. (52) 
Substitution of (52) in (49) yields 



dS OS dK , 1,^ 2/3X. ^ dS 1 

ds dX ds ^ 3 ^ \ dX 2 



1 d^S 



dx^ 

Provided 
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the distinguished hmit of the previous equation for S gives 



A more detailed derivation using a book-keeping small parameter is included in 

Appendix B. Note that the inflection point K{s) is between the position of the 
local maximum of the distribution function for large sizes and the maximum 
size, kM < K. If we change variables from the adaptive time s to the front 
location K, we find an equation in which K scales as X^. Thus it is convenient 
to consider 5" as a function of the new 'time' K and the similarity variable 
^ = X/^/K. Then, 



^ dK 6 ~ 2 ' ^ ^ 



k - K{s) 



(56) 



Equation (55) should be solved with boundary condition S{oo, K) = and an 
appropriate matching condition as ^ — > — oo. 

The solution of (53) is 

K(s) = {?^ + Krf\ (57) 

in which Kq is an arbitrary positive constant to be selected later. The solution 
(28) with a{k) ~ — 3 + 3A;^/^/2 yields the matching condition 

a,^2c(^s + 3-le/') r.2c(3-lK'J'-^K'f') , (58) 



in the overlap region: y/K <^ [K — k) ^ /C, as ^ ^ — cxo. In Appendix C, we 
show that the solution of (55) satisfying boundary and matching condition is 

S{i. K) = ^ / [c{t')r e dt'. (59) 

Clearly, the front solution (59) contributes to the moment M1/3 in (21) for 
small times corresponding to k in the overlap region. Since S is matched to 
a variable outer solution, it is convenient to pick a time tp{t) corresponding 



14 




Fig. 5. Nondimensional size distribution function r/c(t) evaluated using the compos- 
ite solution (61) - (62) (dashed line) and the numerical solution of the full model 
discrete equations (solid line) at the nondimensional times 300, 400 and 500. 

to k in the overlap region to split the time interval (0,t) in (48) into two 
subintervals. For < t' < tp(t), we use the front approximation (50), (52) 
and (59) whereas we use the outer approximation (28) with (45) for times in 
{tp{t),t). The patching time solves 



s{tp)=^p[K{s{t))f' + ^-lK'/', (60) 

for in the overlap region. Since the width of the gaussian in (59) is 76, we 
may choose ^p > Up to the patching time, the contribution of the leading 
edge (59) to the moment M1/3 in (21) is S{^, K) K^/^d^ = 2 K^/^ /q^ [c{t')f dt' . 
Taking into account the approximation (28), we obtain 

tp{t) 

+2[K{sit))Y/^c j [c{t')fdt' = 1. (61) 




Fig. 5 compares the numerical solution of the full discrete model equations to 
the composite solution: 
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= 2k-'^/^c{s{t) - a{k))+e{K - ^p^/K - k) 

mk-K + i,^) e(../3-.V3) dt\ (62) 

V67r(i^V3 _ i^V3) I 

plus (61) for the monomer concentration. We have used the numerical values 
Kq = 0.5 and = a/6 without looking for an optimal fit to the numerical 
solution of the full model equations by varying Kq and C,p- The agreement 
between our theory and the numerical solution of the full model is much better 
than that achieved by using only the outer solution. Although our effective 
small parameter is the reciprocal of the wave front location (and l/X — > as 
t — > cxo), the poor agreement between the outer solution and the numerical 
solution at large sizes precludes a better agreement between the composite 
solution and the numerical solution. 



5 Moment closure 



The composite solution is computationally costly if all we want to calculate 
is the monomer concentration because we need to calculate K{s{t)), tp{t) and 
c(s — a{k)) at each instant of time. A different approximation involving only 
equations that are local in time was used by Schaldach and Wolfer [1] to 
approximate the number densities of monomers, bubbles and helium. Note 
that Equations (21), (23) and (24) would form a closed system of equations if 
we could express M1/3 in terms of Mq and Mi. Physically speaking, M1/3/M0 
and Ml /Mo have the meaning of average bubble radius and average helium 
density. If we impose 



Mo 3 V Afo ; ' KA-kJ 10, \ ) 

we obtain the following closed system of three differential equations [1] 



dMo 
ds 



2c, (64) 



4c+f-^V^'M^/X/^ (65) 



ds \47r 
dc , , , , / 3 \ 

ds 



+ 4c^+(|-) cMl^'M',/'^!, (66) 



plus Eq. (20) relating s and the time. For long times, these equations possess 
the same scaling symmetry as the continuum equations and their solutions 
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Fig. 6. Nondimensional size distribution function rfc(t) of (29) and (45) with 
monomer concentration given by the moment equations (64) - (66) (dashed hne) 
and the numerical solution of the full model discrete equations (solid line). Times 
are as in Fig. 3. 

tend to a similarity solution for long times. A different closure assumption can 
be obtained by combining the idea of preserving the scaling symmetry with a 
closure assumption related to (63), as explained in Appendix D. The resulting 
moment equations have solutions which become self-similar for appropriate 
intervals of time, but these solutions approximate the numerical solution of 
the kinetic equation less well than those of (64) - (66). 

Figure 6 shows a comparison between the size distribution function (29) and 
(45) calculated with the monomer concentration resulting from the moment 
equations (64) - (66) and the numerical solution of the full model. We observe 
that the maximum at large sizes is predicted to occur at much larger sizes. 
Thus the moment equations provide a worse prediction than the outer solution 
(48), and, of course, worse than the composite solution (61) - (62). 



6 Conclusions 



In Section 3, we have found the composite solution (61) - (62) to the discrete 
kinetic equations. The outer approximation is a solution of the continuum 
limit of the discrete equations corrected by the effects of discreteness. The in- 
ner approximation follows from a wave front expansion previously introduced 
for transient homogeneous nucleation [4] . A similar equation for the wave front 
profile was obtained by King and Wattis in the case of the Becker-Doring equa- 
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tions with rate constants having a power law dependence on cluster size [5]. 
While in the case of nucleation the outer approximation was a simple con- 
stant profile in the wake of the wave front, in the present case of growing 
helium bubbles, the outer approximation is a time and size dependent func- 
tion depending on the monomer concentration. The monomer concentration 
satisfies an integrodifferential equation comprising contributions of both the 
outer and the inner solutions. As time increases, the outer solution tries to 
achieve a self-similar form in the variable x = A;/ s^^^ corrected by discreteness 
effects, as suggested by Fig. 2. The inner solution is determined by the posi- 
tion of the wave front K(s{t)) and by the similarity variable ^ = (k — K)/ \fK^ 
which is different from x- Inner and outer solution have an overlap domain 
1 <^ ^ [K — k) <^ K and are patched at a point thereof to yield the 
composite solution. 

We have also compared our approximations to Schaldach and Wolfer's mo- 
ment closure of Section 5 supplemented with the formulas (29) and (45) for 
the size distribution function. Fig. 3 shows that the similarity solution for c(t) 
is consistently above the numerical solution of the model, while the moment 
closure solution is consistently below. Then the maximum of the size distri- 
bution function approximated using the similarity solution (resp. the moment 
closure theory) occurs before (resp. after) the real maximum. In contrast with 
these approximations, the outer solution (48) yields a better match to the 
monomer concentration for large times and to the location of the maximum of 
the size distribution function. Corrected with the boundary layer formula, the 
composite solution (61) - (62) gives the best match to the numerical solution 
of all the approximations described in this paper. 
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A Approximating the sums Ylj=2j 

For < q; < 1 and /c — > oo, these expressions are Riemann sums whose jth 
term equals the area under a rectangle of height (j — 1)^" and basis 1. These 
Riemann sums can be approximated by the integral minus the sum of the 
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areas of triangles of height [{j — 1)" — j"] and basis 1. We thus have 

Ei"" = Ei-" - 1 ^ / ^ + ^ Eb-" - (j - - 1, 

j=2 J=l ^ ^ ,=2 

and therefore 

E^ ~i ^^^^^ — 



^.^2 - a 2(1 -a) 2k'^' 

The function a(/c) is given by the sum with exponent 1/3, which yields (45). 
Similarly, the half- width of the peak is given by (7(/c), 



k 

a{kr ^ f[s- a{k)r Rk{s) ds = R"{0) - [R'{0)f = Er'^' 

J=2 

As explained in the text, the relative width of the peak of Rk{s), aja 
2k-^/'^/^/?, tends to as A; ^ oo. 



B Derivation of the equation for the wave front profile using a 
book-keeping small parameter. 



Let us insert 



(Tfc = S{X\ s*), X* l^lfl j , s* = e'/3^, e ^ 0+, 

in (49) instead of (52). The small parameter e represents location of the wave 
front at a typical large size. We find 



^2/3 



dS ^_,,,dSdK* f{K*y/' £2/3-7 \ 



/ , dS e^T d'S 

X —f^ I u 

\ dx* 2 d{x*f ■ 



from which 



V^s* 3(K*)2/3aX*; 9X* ' ds* 
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^27-1/3 1/3 Q2g 

2 d{x*y ^ 

Provided 

and 7 = 1/2, we find the equation 

dS X" dS {K*f''' d^S 



ds* 3{K*)y^dX* 2 

in the hmit as e — > 0+. The previous two equations are (53) and (54) once we 
revert to the variables s and K. 



C Solution of the equation for the leading front 

Defining J — —dS/d^, we find the following equation for J: 

dJ 1 dj^J) _ 1 d'J 
dK 6 dC 2 9^ ■ 

The Fourier transform of J(^, •) satisfies the hyperbolic equation 

dK 6 dr) 2 ' 

which is readily solved by the method of characteristics in terms of an arbitrary 
initial condition jo(?7o)- Inverting the Fourier transform and going back to the 
function 5", we obtain 



1 

5(C, K) = S{K) + , , , ,, 7^ / -^0(^0) exp 



(e-6 



Q[1-{K,/KY/^] 



in which S{K) and 5*0(^0) are both arbitrary. As ^ — >• —00, the integral can be 
approximated by the Laplace method with the result S{^, K) ~ S{K)-\-Sq{^). 
The matching condition (58) gives Sq ~ 2c(3 - K^'^i - ZKI'^ /2) - S{K), 
thereby yielding 



S(^, K) = ~ / c i^'/'ei + 3 ^ e 6[i-(Ko/K)i/3] 
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Changing variables from to s' — K^^^^i + 3 — 3Kq'^/2, we find 



3if2/3 \ 2 
Kl/6f_3+ff0_+/ 

K) = , " I c{s') e «<-^^^-o^^) d.'. (C.l) 

Tfie ends of the integration interval in this expression are set to and s because 
those are the extremes of the interval over which the monomer concentration 
exists. Changing variables in this formula from s' to the time t', we obtain 
(59). 

The width of the Gaussian in (C.l) is VE. Thus a typical point in the overlap 
region is kp — K — ^p^/K with ^p > -\/6. The corresponding adaptive time is 



and the corresponding time is given by (60). 



D Moment equations following from a closure assumption preserv- 
ing scaling symmetry 



Let us assume that the size distribution function has the scaling form 



^k{s)^ ,77'l f{x), (D.l) 



Mi(s)Zo 

oo 

[ x^r{x) dx. (D.3) 







The definition (22) of the moments together with (D.l) - (D.3) imply 



M^^-^MfMo^-^ (D.4) 

We can use (D.4) for = 1/3 to close the system of equations (21), (23) and 
(24), with the result 
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dMp 

ds 



2c, 



as 



dc 



1/3, .2/3 



c— + 4c' + XcMi'^MQ 
ds 

ds 



dt 
X- 



= c, 



tl/3 



,1/3 ,2/3 • 

n '0 



(D.5) 
(D.6) 
(D.7) 
(D.8) 
(D.9) 



These equations become (64) - (66) if A = (47r/3)-i/3^ To make them com- 
patible with the previously found similarity solution, we consider the reduced 
system given by (D.5) and the approximate equations 



ds c 



(D.IO) 



which have the same scaling symmetry as in the case of the continuum limit. 
Thus they have the similarity solution 



c = C,s-^/\ 



(D.ll) 



for a certain constant Ci. Prom (D.5) and (D.IO) (written as dM^/ds — 1/c), 
we obtain 



Mo^9>CiS^'\ Ml 



whereby (D.IO) yields 



7Ci ' 



(D.12) 



71/4 



(D.13) 



We shall now use Eq. (35) for and (D.l) and (D.ll) - (D.13) to find f. 
Straightforward but lengthy calculations yield 



f{x) — 



I, 



Xm 



/ X 

6xm \xm 

.6/1/3J 



-1/3 



1- — 



X \2/3' 



n-3/4 



XmJ 



(D.14) 
(D.15) 
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With this reduced size distribution function we can check that (D.3) becomes 
Iq — Iq ioT II — whereas it becomes 



7lf loTH/if 



\ 3nl 



12 ' ^^-^^^ 



1/3 



3/v3 km/ AY 

for = 1 and 1/3, respectively. These two last equations and (D.15) imply 
that 



^1/3 = ^^1, XM^l. (D.18) 

Thus (D.18) are required for the reduced size distribution function to be con- 
sistent with the definitions of the 1^. Using (D.16) - (D.18) in (D.9) and (D.13), 
it is possible to show that Ci = -Ro/2 given by (38) and we get the same sim- 
ilarity solution as before. 

Had we used Schaldach and Wolfer's closure A = (47r/3)~^/^, we would have 
obtained the following similarity solution to the reduced system (D.5) and 
(D.IO): 



If we employ time instead of the variable s, (D.19) becomes 



c 



47r\ ^_2/7^_3/7^ ^ / 47r\ ^3/^^^/^^ ^ ^ ^^^^^ 
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